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3D Computer simulations and experiments are employed to study random packings of compress- 
ible spherical grains under external confining stress. Of particular interest is the rigid ball limit, 
which we describe as a continuous transition in which the applied stress vanishes as (0 — (fic) 13 , where 
4> is the (solid phase) volume density. This transition coincides with the onset of shear rigidity. The 
value of 4>c depends, for example, on whether the grains interact via only normal forces (giving rise 
to random close packings) or by a combination of normal and friction generated transverse forces 
(producing random loose packings). In both cases, near the transition, the system's response is 
controlled by localized force chains. As the stress increases, we characterize the system's evolution 
in terms of (1) the participation number, (2) the average force distribution, and (3) visualization 
techniques. 
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Dense packings of spherical particles are an important 
starting point for the study of physical systems as di- 
verse as simple liquids, metallic glasses, colloidal suspen- 
sions, biological systems, and granular matter 0^]. In 
the case of liquids and glasses, finite temperature molec- 
ular dynamics (MD) studies of hard sphere models have 
been particularly important. Here one finds a first or- 
der liquid-solid phase transition as the solid phase vol- 
ume fraction, </>, increases. Above the freezing point, a 
metastable disordered state can persist until <f> — > R cp 
|^[ , where </>rcp is the density of random close packing 
(RCP) — the densest possible random packing of hard 
spheres. 

This Letter is concerned with the non-linear elastic 
properties of granular packings. Unlike glasses and amor- 
phous solids, this is a zero temperature system in which 
the intcrparticle forces are both non- linear, and path (i.e., 
history) dependent. [Because these forces are purely re- 
pulsive, mechanical stability is achieved only by impos- 
ing external stress.] The structure of packing depends 
in detail on the forces acting between the grains during 
rearrangement of grains; indeed, different rearrangement 
protocols can lead to either RCP or random loose packed 
(RLP) systems. 

In the conventional continuum approach to this prob- 
lem, the granular material is treated as an elasto-plastic 
medium [pi. However, this approach has been challenged 
by recent authors |q] who argue that granular packings 
represents a new kind of fragile matter and that more 
exotic methods, e.g., the fixed principal axis ansatz, are 
required to describe their internal stress distributions. 
These new continuum methods are complemented by mi- 
croscopic studies based on either contact dynamics simu- 
lations of rigid spheres or statistical models, such as the 
q-modcl, which makes no attempt to take account of the 
character of the inter-grain forces [^],[| . 

In our view, a proper description of the stress state 
in granular systems must take account of the fact that 



the individual grains arc dcformablc. We report here 
on a 3D study of deformable spheres interacting via 
Hertz-Mindlin contact forces. Our simulations cover four 
decades in the applied pressure and allow us to under- 
stand the regimes in which the different theoretical ap- 
proaches described above are valid. Since the grains in 
our simulations are deformable, the volume fraction can 
be increased above the hard sphere limit and we are able 
to study the approach to the RCP and RLP states from 
this realistic perspective. Within this framework, the 
rigid grain limit is described as a continuous phase tran- 
sition where the order parameter is the applied stress, 
cr, which vanishes continuously as ((f) — (^c) 13 ■ Here cf> c is 
the critical volume density, and (3 is the corresponding 
critical exponent. We emphasize that the fragile state 
corresponding to rigid grains is reached by looking at the 
limit 0^0+ from above. 

Of particular importance is the fact that <j) c depends on 
the type of interaction between the grains. If the grains 
interact via normal forces only fief , they slide and ro- 
tate freely mimicking the rearrangements of grains dur- 
ing shaking in experiments jl]-^. We then obtain the 
RCP value <\> c = 0.634(4) (« RC p). By contrast, if the 
grains interact by combined normal and friction gener- 
ated transverse forces, we get RLP jj] at the critical point 
with (j> c — 0.6284(2) < </>rcp- The power-law exponents 
characterizing the approach to 4> c are not universal and 
depend on the strength of friction generated shear forces. 

Our results indicate that the transitions at both RCP 
or RLP are driven by localized force chains. Near the crit- 
ical density there is a percolative fragile structure which 
we characterize by the participation number (which mea- 
sures localization of force chains), the probability dis- 
tribution of forces, and also by visualization techniques. 
A subset of our results are experimentally verified using 
carbon paper measurements to study force distributions 
in the granular assembly. We also consider in some de- 
tail the relationship between our work and recent exper- 
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iments in 2D Couette geometries [fij 
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FIG. 1. (a) Confining stress and (b) average coordination 
number as a function of volume fraction for friction and fric- 
tionless balls. 

Numerical Simulations: To better understand the be- 
havior of real granular materials, we perform granular 
dynamics simulations of unconsolidated packings of de- 
formable spherical glass beads using the Discrete Element 
Method developed by Cundall and Strack fl^] . Unlike 
previous work on rigid grains, we work with a system 
of deformable elastic grains interacting via normal and 
tangential Hertz-Mindlin forces plus viscous dissipative 
forces . The grains have shear modulus 29 GPa, Pois- 
son's ratio 0.2 and radius 0.1 mm. 

Our simulations employ periodic boundary conditions 
and begin with a gas of 10000 non-overlapping spheres 
located at random positions in a cube 4 mm on a side. 
Generating a mechanically stable packing is not a triv- 
ial task || . At the outset, a series of strain-controlled 
isotropic compressions and expansions are applied until 
a volume fraction slightly below the critical density. At 
this point the system is at zero pressure and zero coordi- 
nation number. We then compress along the z direction, 
until the system equilibrates at a desired vertical stress 
a and a non-zero average coordination number (Z). 

Figure |l|a shows the behavior of the stress as a func- 
tion of the volume fraction. We find that the pressure 
vanishes at a critical 4> c — 0.6284(2). Although we can- 
not rule out a discontinuity in the pressure at </> c — as 
we could expect for a system of hard spheres — our re- 
sults indicates that the transition is continuous and the 
behavior of the pressure can be fitted to a power law form 
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where j3 — 1.6(2). Our 3D results contrast with recent 
experiments of slowly sheared grains in 2D Couette ge- 
ometries 1 11 where a faster than exponential approach 
to <j> c was found, while they agree qualitative with similar 
continuous transition found in compressed emulsions and 
foams |lC]_ 

lb shows the behavior of the mean coordination 
(Z) , as a function of </>. We find 
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where Z c = 4 is a minimal coordination number, and 
9 = 0.29(5) is a critical exponent. At criticality the sys- 
tem is very loose and fragile with a very low coordination 



number. The value of Z c can be understood in term of 
constraint arguments as discussed in |iH ; in the rigid 
ball limit, for a disordered system with both normal and 
transverse forces, we find Z c = D + 1 = 4 Gj]. As we 
compress the system more contacts are created, provid- 
ing more constraints so that the forces become overde- 
termined. 

We notice that <fi c obtained for this system is consid- 
erably lower than the best estimated value at RCP [pj, 
0rcp = 0.6366(4) obtained by Finney using ball bear- 
ings. This latter value is obtained by carefully vibrating 
the system and letting the grains settle into the most 
compact packing. Numerically, this is achieved by al- 
lowing the grains reach the state of mechanical equilib- 
rium interacting only via normal forces. By removing the 
transverse forces, grains can slide freely and find most 
compact packings than with transverse forces. Numeri- 
cally we confirm this by equilibrating the system at zero 
transverse force. The critical packing fraction found in 
this way is <p c = 0.634(4) (« 0r,cp within error bars). 
The stress behaves as in Eq. ([[]) but with a different ex- 
ponent /3 — 2.0(2) (Fig. 0a). At the critical volume frac- 
tion the average coordination number is now Z c = 6 [and 
9 = 0.94(5), Fig. |]b], which again can be understood us- 
ing constraint arguments which would give a minimal co- 
ordination number equal to 2D for frictionless rigid balls 
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FIG. 2. Distribution of forces for different confining 
stresses a obtained (a) in the numerical simulations of fric- 
tion balls, and (b) in the carbon paper experiments. The 
straight solid lines are fittings to exponential forms and the 
dashed lines are fittings to Gaussian forms. In both graphs we 
shift down the distributions corresponding to the two larger 
stresses for clarity, (c) Participation number versus external 
stress for the same system as in (a). 

We conclude that the value 4> c w 0.6288 < </>rcp found 
with transverse forces corresponds to the RLP limit, ex- 
perimentally achieved by pouring balls in a container but 
without allowing for further rearrangements [Q. Exper- 
imentally, stable loose packings with tj> as low as 0.60 
have been found Q. In our simulations, C lower than 
0.6288 can be obtained by increasing the strength of the 
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tangential forces. This is in agreement with experiments 
of Scott and Kilgour [fisj who found that the maximum 
packing density of spheres decreases with the surface 
roughness (friction) of the balls. 

While previous studies characterized RCP's and RLP's 
by using radial distribution functions and Voronoi con- 
structions ||], we take a different approach which allow 
us to compare our results directly with recent work on 
force transmissions in granular matter. Previous studies 
of granular media indicate that, for forces greater than 
the average value, the distribution of inter-grain contact 
forces is exponential (^J^]. In addition, photo-elastic vi- 
sualization experiments and simulations Jl6|,p| show that 
contact forces are strongly localized along "force chains" 
which carry most of the applied stress. The existence 
of force chains and exponential force distributions are 
thought to be intimately related. 

Here we analyze this scenario in the entire range of 
pressures: from the <f> c limit and up. Figure |^a shows 
the force distribution obtained in the simulations with 
friction balls. At low stress, the distribution is exponen- 
tial in agreement with previous experiments and models. 
When the system is compressed further, we find a gradual 
transition to a Gaussian force distribution. We observe 
a similar transition in our simulations involving friction- 
less grains under isotropic compression. This suggests 
that our results are generic, and do not depend, qualita- 
tively, on the preparation history or on the existence of 
friction generated transverse forces between the grains. 

Physically, we find that the transition from Gaussian 
to exponential force distribution is driven by the local- 
ization of force chains as the applied stress is decreased. 
In granular materials, with particles of similar size, local- 
ization is induced by the disorder of the packing arrange- 
ment. To quantify the degree of localization, we consider 
the participation number II: 

ri=f M X>/) . (3) 

Here M is the number of contacts between the spheres, 
(Z) = 2M /N is the average coordination number, and N 

is the number of spheres. % = fi/Ylj=i fjt wnere fi is 
the magnitude of the total force at every contact. From 
the definition (SI), H = 1 indicates a limiting state with a 
spatially homogeneous force distribution (<^ = 1/M, Vi). 
On the other hand, in the limit of complete localization, 
n w 1/M -> and M -» oo. 

Figure |^c shows our results for n versus a. Clearly, 
the system is more localized at low stress than at high 
stress. Initially, the growth of n is logarithmic, indicat- 
ing a smooth derealization transition. This behavior is 
seen up to a w 2.1 MPa, after which the participation 
number saturates to a higher value: 

n(cr) oc log(er) [a <2.1 MPa] , , 

n(cr) w 0.62 [cr >2.1 MPa] ( > 



This behavior suggests that, near the critical density, the 
forces are localized in force chains sparsely distributed in 
space. As the applied stress is increased, the force chains 
become more dense, and are thus distributed more ho- 
mogeneously. 

How might we expect the participation number to de- 
pend upon other system parameters when the forces are 
transmitted principally by force chains? In an idealized 
situation, the system has Npc force chains, each of which 
has N z spheres. Each sphere in a force chain has two 
major load bearing contacts, which loads must be ap- 
proximately equal. In the lateral directions, roughly four 
weak contacts are required for stability. These contacts 
carry a fraction a < 1 of the major vertical load. All 
other contacts have fi » 0. Under these assumptions, 

2 (1 + 2a) 2 N FC N Z 2 (1 + 2a) 2 
(Z) (1 + 2a 2 ) N ~ (Z) (1 + 2a 2 ) ' 1 ' 

The last inequality becomes an equality iff all the balls 
are in force chains. From our simulations at large pres- 
sure a w 2/5, so at (Z) « 8, n « 0.62, which implies that 
the system has been completely homogenized. Although 
Eq. <S) is oversimplified, we believe that the change in 
slope in Fig. ||c is emblematic of the complete disappear- 
ance of well-separated chains. 
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FIG. 3. Example of percolating force chains for the same 
system as in Fig. |^a: (a) near <j> c and (b) away from (f> c at 
large confining stress. The color code of the chains is accord- 
ing to the total force in N carried by the chains. 

The localization transition can be understood by 
studying the behavior of the forces during the loading 
of the sample. Clearly, visualizing forces in 3D systems 
is a complicated task. In order to exhibit the rigid struc- 
ture from the system we visually examine all the forces 
larger than the average force; these carry most of the 
stress of the system. The forces smaller than the average 
are thought to act as an interstitial subset of the system 
providing stability to the buckling of force chains . 
We look for force chains by starting from a sphere at the 
top of the system, and following the path of maximum 
contact force at every grain. We look only for the paths 
which percolate, i.e., stress paths spanning the sample 
from the top to the bottom. In Fig. || we show the evo- 
lution of the force chains thus obtained for two extreme 
cases of confining stress. We clearly see localization at 
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low confining stress: the force-bearing network is con- 
centrated in a few percolating chains. At this point the 
grains are weakly deformed but still well connected. We 
expect a broad force distribution, as found in this and 
previous studies. As we compress further, new contacts 
are created and the density of force chains increases. This 
in turn gives rise to a more homogeneous spatial distri- 
bution of forces, which is consistent with the crossover to 
a narrow Gaussian distribution. 

Experiments: Some of the predictions of our numer- 
ical study can be tested using standard carbon paper 
experiments || , which have been used successfully in the 
past to study the force fluctuations in granular packings. 
A Plexiglas cylinder, 5 cm diameter and varying height 
(from 3 cm to 5 cm), is filled with spherical glass beads 
of diameter 0.8 ± 0.05 mm. At the bottom of the con- 
tainer we place a carbon paper with white glossy paper 
underneath. We close the system with two pistons and 
we allow the top piston to slide freely in the vertical direc- 
tion, while the bottom piston is held fixed to the cylinder. 
The system is compressed in the vertical direction with 
an Inktron™ press and the beads at the bottom of the 
cylinder left marks on the glossy paper. We digitize this 
pattern and calculate the "darkness" ||] of every mark 
on the paper. To calibrate the relationship between the 
marks and the force, a known weight is placed on top 
of a single bead; for the forces of interest in this study 
(i.e., from w 0.05N to 6 N), there is a roughly linear re- 
lation between the darkness of the dot and the force on 
the bead. 

We perform the experiment for different external 
forces, ranging from 2000 N to 9000 N, and different 
cylinder heights. The corresponding vertical stress, a, 
at the bottom of the cylinder ranges between 100 KPa 
and 2.3 MPa (as measured from the darkness of the dots). 
The results of four different measurements are shown in 
Fig. ||b. For a smaller than « 750 KPa, the distribution 
of forces, /, at the bottom piston decays exponentially: 

P(f) = if)' 1 cxp [-//</>], [ <r < 750 KPa], (6) 

where (/} is the average force. When the stress is in- 
creased above 750 KPa there is a gradual crossover to a 
Gaussian force distribution as we find in the simulations. 
For example, at 2.3 MPa we have 



P(/)cxexp -k 2 (f-f o y 



[cr=2.3 MPa]. (7) 



where kf a w 1, and therefore (/) « f . Similar results 
have been found in 2D geometries pM . 

Discussion: In summary, using both numerical simu- 
lations and experiments, we have studied unconsolidated 
compressible granular media in a range of pressures span- 
ning almost four decades. In the limit of weak compres- 
sion, the stress vanishes continuously as (</> — 4> c )^ , where 
(f> c corresponds to RLP or RCP according to the exis- 
tence or not of transverse forces between the grains, re- 
spectively. At criticality, the coordination number ap- 
proaches a minimal value Z c (=4 for friction and 6 for 



frictionless grains) also as a power law. Our result Z c = 6 
agrees with experimental analysis of Bernal packings for 
close contacts between spheres fixed by means of wax , 
and our own analysis of the Finney packings Q using 
the actual sphere center coordinates of 8000 steel balls. 
However, no similar experimental study exists for RLP 
which could be able to confirm Z c = 4. A critical slowing 
down — the time to equilibrate the system increases near 
4> c — and the emergency of shear rigidity (to be discussed 
elsewhere) is also found at criticality. The distribution 
of forces is found to decay exponentially. The system 
is dominated by a fragile network of relatively few force 
chains which span the system. 

When the stress is increased away from 4> c to the point 
that the number of contacts has significantly increased 
from its initial value Z c we find: (1) the distribution of 
forces crosses over to a Gaussian (2) the participation 
number increases, and then abruptly saturates and (3) 
the density of force chains increases to the point where 
it no longer makes sense to describe the system in those 
terms. Our simulations indicate that the crossover is 
associated with a loss of localization and the ensuing ho- 
mogenization of the force-bearing stress paths. 
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